%% xiaoyeyimier
% Date : 2021-10-12
% Version : 1.0.2
% Fuction : Receive GPS data via BladeRF
% Contact : xiaoyeyimier(Wechat)
% Instructions :
%
%   The number of reads in memory at a time is an issue to consider. It is 
%   set to read memory data every 3 seconds and store it in a file currently.
%-------------------------------------------------------------------------%
%% Date : 2022-09-24
% Version : 1.0.1
% Try adding some ionospheric scintillation measurements Based on BladeRF_GPSL1 1.0.2
%-------------------------------------------------------------------------%
close all;
clearvars all;
clc;addpath('CN0Estimate');
%% Initial
%-------------------------------------------------------------------------%
[DevSettings,LogID] = SetBladeRF();                                        % set hardware
[AcqSettings,TrackingSettings] =preRunStandardReceiver(DevSettings);               % set software
addpath(DevSettings.LibPath);                                              % all lib to current path
%=========================================================================%
%% Start BladeRF 
BladerfFindFlag = true;
if DevSettings.RealTimeProcess
    try
        Dev = bladeRF();                                                   % Try linking to BladeRF
    catch
        fprintf(LogID,'\r\n%s!!\r\n',upper('BladeRF failed to link'));
        warning('BladeRF failed to link,Using a historic data instead!');
        BladerfFindFlag = false;% If the current device is offline, the device connection will be skipped
    end
    if BladerfFindFlag
        fprintf(LogID,'%s!!\r\n',upper('BladeRF successfully linked'));
        try
            [NameI,NameQ] = GPS_bladeRF(Dev,DevSettings);
            clear Dev;
            delete(gcp('nocreate'));
        catch
            fprintf(LogID,'\r\n%s!!\r\n',upper('Data reading failed'));
            error('Data reading failed')
        end
        fprintf(LogID,'%s!!\r\n',upper('Data reading success'));
        fclose(LogID);
    end
end
%% Read Data 
%-------------------------------------------------------------------------%
Times2Read = 15; % Times to read [ms] (for debug and acq)
Samples2Read = DevSettings.Fs*Times2Read*1e-3; % cal samples 2 read
if ~DevSettings.RealTimeProcess || ~BladerfFindFlag % read a realtime data
    FidI = fopen([DevSettings.HistoryDataPath,'\',...
                  DevSettings.HistoryDataName,'Idata.',...
                  DevSettings.FileType],'r');
    FidQ = fopen([DevSettings.HistoryDataPath,'\'...
                  DevSettings.HistoryDataName,'Qdata.',...
                  DevSettings.FileType],'r');
else                                                                       % read a history data
    FidI = fopen(NameI,'r');
    FidQ = fopen(NameQ,'r');
end
fseek(FidI,Samples2Read*16*0.5,'bof');                                      % skip some unstable data from the beginning
fseek(FidQ,Samples2Read*16*0.5,'bof');                                      % skip some unstable data from the beginning
dataI = fread(FidI,[1 Samples2Read],DevSettings.dataType);                 % read I data
dataQ = fread(FidQ,[1 Samples2Read],DevSettings.dataType);                 % read Q data
longSignal = dataI+1j*dataQ;                                               % Combine IQ to a complex signal

%% pulse blanking
%-------------------------------------------------------------------------%
if TrackingSettings.IsBlanking
    Ex = mean(abs(longSignal)); % Compute the expectation of the Rayleigh distribution
    th = sqrt(2/pi)*sqrt(-2*log(TrackingSettings.Pa))*Ex; % Calculate thresholds from expectations
    res = longSignal(abs(longSignal)>th);
    longSignal(abs(longSignal)>th)=th*res./abs(res); % pulse blanking
end
% if TrackingSettings.IsBlanking
%     longSignal=longSignal./abs(longSignal);                                      % pulse blanking
% end
%% plzot time domain figure
%-------------------------------------------------------------------------%
dataI = real(longSignal);
dataQ = imag(longSignal);
t = Times2Read*0.5+(1:Samples2Read)/DevSettings.Fs*1e3;  
subplot(2,1,1);
plot(t,dataI,'r');
xlabel('Time[ms]')
subplot(2,1,2);
xlabel('Time[ms]');
plot(t,dataQ,'b');

%% Acqusition
%-------------------------------------------------------------------------%
acqResults = acquisition_bladerf(longSignal,AcqSettings);
if (any(acqResults.carrFreq))
        channel = preRun(acqResults, AcqSettings);
        showchannel(channel, AcqSettings);
else
    % No satellites to track, exit
    disp('No GNSS signals detected, signal processing finished.');
    trackResults = [];
    return;
end

%% Track the signal
%-------------------------------------------------------------------------%
startTime = now;
disp (['Tracking started at ', datestr(startTime)]);
[trackResults, channel] = tracking_BladeRF(FidI,FidQ,channel,TrackingSettings); 
[trackResults] = getS4(trackResults,channel);
showTrackResults(trackResults,channel);
fclose all;